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Path-integral approach to the tight-binding polaron is extended to multiple optical phonon modes 
. . . of arbitrary dispersion and polarization. The non-linear lattice effects are neglected. Only one 

1/^ ' electron band is considered. The electron-phonon interaction is of the density-displacement type, 

' but can be of arbitrary spatial range and shape. Feynman's analytical integration of ion trajectories 

^— »' , is performed by transforming the electron-ion forces to the basis in which the phonon dynamical 

matrix is diagonal. The resulting polaron action is derived for the periodic and shifted boundary 
conditions in imaginary time. The former can be used for calculating polaron thermodynamics while 
the latter for the polaron mass and spectrum. The developed formalism is the analytical basis for 
numerical analysis of such models by path-integral Monte Carlo methods. 
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In the last two decades, quantum Monte Carlo (MC) simulation methods proved to be a powerful theoretical tool in 
analyzing models with strong electron-phonon interactions. The pioneering works by Hirsch and Fradkin on the many- 
polaron one-dimensional Holstein modeli and the Su-Schrieffer-Heeger polaron- , by de Raedt and Lagendijk on the 
^ Holstein polaron'^ and bipolaron^, and by Alexandrou and Rosenfelder on the Frohlich polaronA, were followed more 
I _ recently by the continuous-time path-integral MC^, diagrammatic MC^'^, and the Lang-Firsov MC^ methods. The 
' O ' new methods have expanded the list of calculable polaron properties, both static and dynamic. However, expansion in 
^ I terms of analysis of more complex polaron models has been slow. The bulk of new literature is still devoted to the two 
^ , canonical polaron models: Frohlich model of the ionic solid^'^ and Holstein model of molecular crystalii. While these 
models are important, the new MC methods are powerful and versatile enough to study others. Notable exceptions 
from this trend are the investigations of the Su-Schrieffer-Heeger model^, Jahn- Teller polaron^^, and of the tJ model 
with electron-phonon interactioni^. 

We recently started a series of generalizations of the continuous-time path integral MC method away from the 
standard Holstein interaction with localized phonons on cubic lattices. A long-range electron-phonon interaction was 
QQ investigated in Refs. ^md different Bravais lattices in Ref. In this paper, we take the first step toward 

; generalizing the method on dispersive and multimode phonon systems. A critical element of the method is analytical 
^-H ' integration over ion trajectories in the path-integral expression for the polaron partition function. This leads to a 
, polaron action, expressed as a functional over the electron path, which includes all the effects of the phonon dynamics 
■ and electron-phonon interaction. For the Einstein phonons, the integration results in the single-oscillator formula 
derived by Feynman-^'^^. In more complex cases, the electron-ion forces must be first transformed to the basis of the 
lattice normal modes, then the normal modes integrated out using the Feynman formula, and then the final result 
' transformed back to the orignal basis. To our knowledge, the only investigation of this kind in the polaron MC 
literature was done by de Raedt and Lagendijk^^ who studied the one-dimensional Holstein polaron with dispersive 
(~| optical phonons by the discrete-time path- integral MC. They found that the transition into a self-trapped state always 
O , takes place for for dispersive phonons as well as for dispersionless phonons but the transition point depends on the 
O ' details of the phonon dispersion. Their derivation now needs to be repeated for continuous imaginary time and 
^ ' generalized to multiple phonon modes. Recently, ZoliiS, calculated the phonon path integral for the Holstein model 
in the semiclassical approximation where the retarded nature of the phonon action was neglected. 
. _ Phonon dispersion is expected to affect polaron properties, most notably the polaron effective mass and hence 
^ . mobility. The forces acting between the ions spread out the deformation even if the electron-phonon interaction is 
■ " ' local. When the electron hops between the lattice cells, partial deformation already exists, which increases the phonon 
overlap integral. Thus in general, dispersion reduces the polaron mass. In this sense, the dispersion is analogous to a 
long-range electron-phonon interactioni^. 

The purpose of the present paper is to carry out general phonon integration to enable MC studies of novel polaron 
models. The main analytical result is given below in equations (|34|l - H36f) and (|31|) - (|32|l . The expressions are rather 
complex, however, and the practical impact of the formalism will depend on the efficiency of numerical methods used 
to evaluate the polaron action at each MC update. This is likely to become the subject of future work. 
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II. SHIFT PARTITION FUNCTION 



Consider the general problem of a lattice electron interacting with an arbitrary system of ions via a density- 
displacement interaction. In the path integral formalism.^^, the electron is described by an imaginary-time path r(T), 
where r is the discrete lattice coordinate at "time" r, and < t < (3 = {kBT)~^. In this paper, only one Wannier state 
per init cell is assumed. The ion system is represented by p > 1 continuous degrees of freedom (displacements) 
where n numbers the unit cells, and 1 < s < p. The mass of s-th degree of freedom is m^. Ions interact with each other 
via pairwise force constants Wns,n's', which are functions of the difference n — n' only. Cubic and higher terms in the 
lattice potential energy are neglected. A displacement £_ns is acted upon by the electron with force — fns[Tc{T)], 

so that the interaction potential is — /ns(''')^ns- The imaginary time action of such a system is given by 

A[v{T),^Ur)] = r[r(T)]+ VK„,(r)]+Ae-ph[r(r),ens(T)] 

= r[r(r)]-^r^^^4l^dr-i ^ / W.'Cn.MCn's'W dr -I- / /„,[r(r)]U(r)dr. (1) 



Here the first term T represents the kinetic energy of the electron. It is not the subject of the present work, see 
instead Refs. 0^3- Note that the action A is dimensionless. The statistical weight of a configuration is given by 
exp {A) . The central object of the polaron path- integral quantum Monte Carlo method is the shift partition function 

^Ar = E / ITdensO / I?r(T) (r) e^^''^^^ ^ ^^'^^^^ . (2) 

ro J ns -'r(0)=ro;5„,(0)=C„,o 

Here the inner path integral is taken under shifted boundary conditions in imaginary time. This means that the entire 
real-space configuration at the end point t — (3 (which comprises the electron and phonon coordinates) is identical 
to the configuration at the end point t — Q up to a parallel shift by a real space vector Ar. The outer ordinary but 
multidimensional integral is over all possible configurations at r = 0. 

Since the ion coordinates enter the action linearly and quadratically, path integration over ij) in can be done 
analytically. We first calculate a more general path integral 

Here, in contrast to eq. the end coordinates of ions are not restricted in any way. is a function of the end 
coordinates ^nsOj ^ns/3j and a functional of the whole electron path. Path integration is performed by a method 
devised by Feynmanii. Variation of the exponent with respect to ^ns{T) subjected to fixed boundary conditions 
yields a Euler-Lagrange equation for the stationary path (r) 

^ LW -Y.^ns,n's'U's'{T) + /n,(T) = . (4) 
i\' s' 

This is a system of pN coupled second order ordinary inhomogeneous differential equations for functions ^(r) with 
boundary conditions 

?"ns(0) =^nsO, (5) 

Us{P)=Usp. (6) 

The inhomogeneous terms fns{T) are functionals of the electron path v{t) and as such are arbitrary functions of 
imaginary time r. When the solution of eq. Q is found, W is calculated by performing a path shift in ||2Jl: 

Us{t) = Us{t) + Vr^sir) , (7) 

VnsiO) = Vnsm=0. (8) 

Upon substitution of {T)) in Q the action Aph + ^c-ph separates in three terms: Ai that depends only on ^, A2 that 
depends only on 77, and the mixed term A3 that depends on both ^ and rj. The mixed term vanishes identically by 
virtue of the classical equations l@J. W becomes 



/•'7n.(/3)=0 
Jr)„s(0)=0 



3 



The path integral over 77 (r) does not depend on any dynamical variables and therefore can be assigned a multiplicative 
constant Z^^^. It represents the incomplete thermodynamic partition function of freely vibrating ions. [The partition 
function will become complete after the final integration over the end ion coordinates in The functional Ai in 
eq. jnjl has the following form: 



2r 



e„.(0)e„s(0)-e„,(/?)C„.(/3) +-/ Us{T'}LAr')dr 



(10) 



As soon as the solution of system (0J is known, W is calculated explicitly using eqs. © and (jTUIl . 



III. SOLUTION OF CLASSICAL EQUATIONS 

In this section we solve the classical equations under the boundary conditions |(SJ|-©. The strategy is to 
transform to normal coordinates, in which the system is diagonal, and then apply the formula for a single 
harmonic oscillator. The first step is a Fourier-mass transformation 

en.(r) = -}^-=^e , (11) 
ak.(r) = V^^ens(T)e-'''", (12) 

n 

where TV is the number of unit cells and k spans the Brillouin zone of the reciprocal lattice. Substituting into 
multiplying by e"^^ " and summing over n, equation Q is transformed into the following: 

-ak.(r) +;i2^w;,v(k)ak.'(T) = ^Ee""'"/-^ ' (13) 
^.v(k) = <,,(k) = ^^'^"""'^ • (14) 

n — n' ^ 

The new p x p matrix w(k) is Hermitian. The latter property is a consequence of the factor in the above 

transformations. This is of course a standard device of the classical theory of lattice vibrations^". The dimensionality 
of li'(k) is frequency squared. On the next step, consider an eigenvalue problem for the matrix K;(k): 

Ywss'{k)ul^ = ujl^u^^. (15) 

s' 

Since w is Hermitian, the eigenvalues are real. We also assume that uj'^^ are positive which implies no structural 
phase transition. Further, since w*^, (— k) = Wss'Ck-) then w_kA = ^kx- The index A numbers different eigenvalues, 
that is phonon modes. The eigenvector u^^ defines the phonon polarization. The eigenvectors are orthogonal, 

,k*„,k 



Ss''^sa'"sA' ^ ^>^>^' ■ k-points of high symmetry, the eigenvectors can be made orthogonal by a proper procedure. 

"'sA — "sA- 



In addition, u^^* — u^^. Next, expand the functions ans(r) in the basis of u: 



aksir) = Ckxir) u^^ , (16) 

A=l 
P 

CkA(T) =Eak.(T)u^^ (17) 



where CkA(T) are new complex functions of imaginary time. Substituting (|lt)f) into (|13|l . using the definition H15|l . 
multiplying the resulting equation by u^y and summing over s one obtains 



CkA(r) + (?i^kA)'ckA(r) = V ^ V e-^'^"/ns(r) . (18) 
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Thus, the system of equations is fully diagonalized, with a transformed inhomogeneous term on the right hand side. 
The problem is reduced to the well studied case of one harmonic oscillator under an arbitrary external force. Before 
applying the corresponding formulas the boundary conditions must be recalculated in the new basis. Combining 
eqs. (O, 10, (|12|l . and (|17|l . the transformed boundary conditions read 

CkA(O) - V^"^A Uso e-'"" , (19) 

ns 



(20) 



To proceed further we recall the following result from the harmonic oscillator theory. A real function H{t) satisfying 
the equation 



-H{t) + {hLUuxfH{T) = F{t) , 
with the boundary conditions H{Q) = and H{P) = Hp, is given by the expression 

sinh?icJkA(/3- t) sinh^WkAT , f'^ r< i '\t?i '\a ' 
H(t) = Ho —— h Hp . + / GkA(T, T )F(t ) dr 



sinh fiuj-k_\j3 sinh fiuj-k.xfi 



where 



sinh huj\^\T ■ sinh huJk\{/3 — r') ; t < t' 



kA(,T, T - ^.^j^ hwkxP \ sinh huji^x{f3 - r) • sinh Tiuj^^xt' ; r > r' 



(21) 



(22) 



(23) 



is the Green's function of the equation (|22|1 under the zero boundary conditions. Separating eqs. H18l) - (|20ll into real 
and imaginary parts and applying formula 122f) one obtains 

, , sinh?i(jJkA(/3 — ■r) , w* -,kn sinh?ia;kAT , k* -,kn 



sinh fkokx/S 



dr'GkA(r, r') ^ ^ ^ e-''"/n.(r') . 



(24) 



Restoration of the classical displacement ^„s from eq. (|ll|l yields the complete solution of the original Cauchy problem 



- 1 sinh hujkx{p ~ 

~ iV^ sinhncc-kA/3 ^ 

kA n's' 

1 sinh huj\^\T 



^ i- Blllll/H^kA' / '"-a' 

AT sinhnwkA/3 ^,\ms 
+ ^Erdr'GkA(r,r')E 



,,k* „,k ik(n-n') 
"s'A^sA ^ 



^msTUs' 



k* k ik(n— n'} 



/nv(r'). 



(25) 



Note that because 



-k* -k -ik(n-n') _ „,k k* -j:k(n-n') 



".'A ".A e 



7/''* p*k(n-n') 



(26) 



the summation over k in equation H25() makes the displacements purely real, as expected. Finally, substituting the 
classical solution in eq. I|1U|) . after some algebra, one obtains the polaron action Ai as a function of the end ion 
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coordinates and a functional of the forces 



^l[CnsO,^ns/3; /ns (''')] 



1 



27V h sinh Tiuj-\^\I3 



y — 



ms'ms 



3? 



k* k zk(n— n') 



X 



nsn s 

X {— (CnsoCn's'O + ^ns/sf n's'/s) COSh TlWkA/? + CnsoCn's'/3 + Cns/3'?n's'o} 

dr 



kA nsn's' 







kA nsn's' 



sinh huJiixP 



, sinh hujy^xT' 



-£.nsafn's'{T') + \ ^^n's'o/ns (t') 



dr' ^ 
Q sinh huj]ixl3 



-^nspfn's'ir') + W ^-Cn's'/3/ns(T') 



1 X - X - h 

kA nsn's' ^ 



7/"^ p«k(n-n') 



/3 /-/J 



dTdr'GuA(T,r')/ns(T)/n's'(T'). 



(27) 



^0 



It should be noted that the real part symbol 3? is optional here, the reality of the action is guarantied by summation 
over k and the above-mentioned symmetry properties of the integrand. Equations © and (|27|l provide an explicit 
solution for the quantity W defined in eq. Q. 



IV. 



POLARON ACTION 



In accordance with the general ideas of the continuous time polaron Monte Carlo, the end points in the action H27|) 
must be related by = Cn-Arso- Then W — Z'^t^ ■ e^^ should be integrated over all ^nso from minus to plus infinity, 
see eqs. 0, and (Q. To perform the integration, it is convenient to employ the same series of transformations 
that enabled calculation of the path integral, that is 



"SnsO — -77 



kA ^ 



:U^Ae'''"gkA-^ 



kA ^ 



,,k*„-ikn * 
:"sAC 5kA 



(28) 



Cns/3 - .fP, 



^«^Ae*("-^"-^5kA ^ 
iV ^ — ' \/ms iV — ' -v/ms 

kA V * kA 



-E — 

N ^ . frn 



^k*g-.k(n-Ar)^*^ 



(29) 



Here ^kA = &kA + *rfkA is a new complex variable. Substituting the expansions into eq. (I27|l one obtains the action as 
a function of amplitudes 6 and d 



^i[&kA,dkA;/ns(r)] = ^E 



t^kA 



N ^-^ h sinh hujit.xl3 

kA 



(cosh;iWkA/3 - coskAr) {-b^x - dix) 



^ E ^-^A E 7^ [-'Ae^*-"] SkAns + 5? 
kA ns V ^ 

E ^kA E 7= [""a^*"] SkAns + 5 
kA ns ^ 



2N ^-^ ^-^ ./m,^m 



^s'A"sA e 



kA nsn's' 



^k^gik(n-A, 

Jo 



where 



B 



kAns 







sinh hujk\{f3 - t') ^ 

sinh hiUkxP ' 



_ /"^^^/ Sinh^tcJkAT^ ^ , ,A 

L-kAns = y fc, ,. , Q J^^y^ ) 



sinh huj]^xP ' 



CkAnsj 
CkAns I 

dTdT'GkA(r,r')/„s(T)/„-,-(T') , (30) 

(31) 
(32) 



and 3? and 3 denote the real and imaginary parts, respectively. The final step is gaussian integration of over 
the real b and d from minus to plus infinity. The reality of ^ and the property u7^* = u^^ ensure that 6-kA = ^kA 
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and d_kA = — c^kA- Thus the domain of integration is hmited to half of the Brillouin zone. The integration is 
straightforward, and leads to the result const ■ e^^' where const completes the partition function of free phonons 
while the polaron action is given by 



AAr[fns{r)] 



N ^ 



h sinh huJif^xP 



E 



1 



4ajkA (cosh hcui^xd — cos kAr) ^-^ ■y./rrisms 

kA nsn's' ^ 



, k ?kn k* — zkn' 



+ —V V ^ 

9/V ^ /-^ 



3? 



k ikn k* — ik(n' — Ar) 



B 



kAnsC'kA n 



kAns^kAn's' 



3? 



ik(r 



-ik(n'-Ar) 



CkAnsCkA: 



n' s' J' 



kA nsn's 



-f Jnisnis' 

/ c/ V 



l/"^* 11^ pjk(n-n') 



/3 /-/J 



dTdr'GkA(T,T')/ns(r)/n's'(T') 



(33) 



^0 



To simplify this expression consider first the case Ar = 0. This corresponds of course to the periodic boundary 
conditions in imaginary time. Taking into account the explicit forms of B, C, and G, the periodic action can be 
brought into the form 



^o[/ns(r)] = l^ ^ 



N ^ 4cjkA Jo Jo 



^rP^ coshac.kA(|-|r-T'|) ^ 
drdr ^ > 

nsn's 



1 



sinh ?ia;kA 



iK p^k(n-n') 



'■sA^s'A': 



/ns(T)/n's'(T') 



(34) 

Comparing this result with an analogous expression for independent oscillators^ shows that eq. (|34ll involves the 
forces transformed to the basis in which the phonon subsystems is diagonal. This is of course exactly what is expected 
intuitively. 

Returning to the general action (|33|l . one can calculate the correction to the periodic action by taking the difference 
between Aat and Aq. The general expression is rather complicated. However, nonzero Ar are needed only in 
calculation of the polaron mass and spectrum, which are the properties of zero temperature partition functions. Thus 
the action correction is needed only in the /3 ^ oo limit. Then cos kAr is negligible compared to cosh?iwkA/5, and 
tanh TiLj\^\(3 ~ 1 with exponential accuracy. Thus the correction becomes 



AylAr(/3^ OO) 



1 X ^ h X ^ 

2^ zL,,„_, 2^ 



1 



N 



4cJkA , , Jnisins 

kA nsn's' ^ 

1 



"sA"s'A 



gik(n— n'+Ar) _ gik(n- 



-Ar) 



ik(n-n'; 



1 X - fl X - 

9,,,,., 



N ^ 2wkA ^ , Jm~m^ 

kA nsn s ^ 



''sA^s 



'A \^ 



GkAns^kAn's'l = 

SkAnsGkA 



+Ar) _ ^ik{n-n') 



(35) 



To summarize the results, the shift partition function |(5J) reduces to a path integral over only the electron coordinates 



^Ar = ^phE / ^^(^) e^[--(^)l+^°["-(^)l+^^^^["-(^)l . (36) 

r„ ^r(0)=ro 

Here = HkAp ^i-'^^(5^/''^kA)]~^ is the full partition of free phonons. It is a multiplicative constant that does not 
influence polaron dynamics and drops out of the configuration weight in Monte Carlo simulations. T is the action 
piece originating from the kinetic energy of the particle. On the lattice and in continuous imaginary time treatment 
of the kinetic energy is non-trivial as has been described elsewhere^'. The periodic action is given by the formula 
()34|1 which is valid for any temperature. This action can be used to calculate equilibrium polaron thermodynamics: 
free energy, specific heat, static correlation functions, and so on. The action correction IS.A is due to shifted boundary 
conditions (non-zero Ar) and is given by equations (|35|l and (|31|I - H32(I . The formulae are valid in the low temperature 
limit, i.e. the difference between eq. ()35|l and the exact formula is 0{e^^'^'^^^). Equation H35() can be used to calculate 
the polaron mass, spectrum, and density of states. If for some reason the full Ar ^ action is needed at a finite 
temperature, then the exact expression (|33|l should be used. 



V. SPECIAL CASES 



The formulae derived in the preceeding section are general but can be difficult to evaluate. In this section, simplified 
expressions for several special cases are obtained that can be used in practical calculations. 
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A. Independent oscillators 

The independence of oscillators implies that the force constants are diagonal in the unit cell index and in the degree 
of freedom index: Uns,n's' = ^sS^n'^ss' ■ Then the matrix elements ^^^'(k) = {vs/ms)6ss' are real and independent 
of the momentum, see eq. (|14ll . The eigenvalue problem H15|l becomes diagonal with eigenvectors Us\ = 6s\ and 
eigenvalues w| = v\/m\, where 1 < X < p. Neither u nor uj depend on k. Substitution into the periodic action H34() 
leads to the following simplifications: (i) Summation over k yields NSnn'] (ii) The emerged delta function removes 
the summation over n'; (iii) Summations over s and s' are performed using the explicit form of the eigenvectors u; 
(iv) The summation index A is changed back to index s to unify the notation. As a result, the periodic function takes 
the form: 

MU.ir)i=±j:^ /rd.d.. -""-'^- I; (37) 

t=i n ^^^n^s Jo Jo smh hcjs | 

where uj^ — yj'vjjrnl. After the same transformations, the action correction (|35f) becomes 

p * 

AA[/„,(t)] =Y.Y1 (Cn+Ar,. - C„.) , (38) 

ZiVgTflg 

S—1 11 

C„.^ / dr' \ Us{t'). (40) 

Jo smhAWsP 

These expressions can be used in analysis of isotropic and anisotropic spherical oscillators in two, three, and higher 
dimensions. 

B. One phonon branch 

Substantial simplifications take place in the case of only one phonon degree of freedom per unit cell with mass m 
and force function w(n — n'). It should be emphasized that this does not restrict the dimensionality of the lattice in 
any way. Rather, the electron interacts predominantly with one type of ion displacement. Mathematically, the indices 
s and A accept only one value 1, and therefore can be omitted hereafter. The eigenvalue matrix equation (|15|1 reduces 
to one linear equation with the solution = \ and 

a;2=u;(k) = lyw(n)e'''". (41) 

TO 

n 

Thus, the phonon frequency is explicitly known. Note that cj^ is automatically real because w(— n) = f (n). These 
results transform the action as follows 



TV ^ 47najk Jo Jo sinh Uu^^^ 



(42) 



AA[/„,(t)] = ^ E E {cos[k(n - n' + Ar)] - cos[k(n - n')]} SknCkn' , (43) 



k 



where B and C are given by H31(l - (|32|l with indices s and A removed. The derived action should be used, e.g., in 
analysis of the polaron on a linear chain of one-dimensional coupled oscillatorsi^. 
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C. Holstein interaction 



The Holstein local interaction is one of the most popular polaron modelsii. At the same time, it is usually studies 
in combination with an Einstein set of decoupled one-dimensional oscillators. Here we will derive the action for the 
Holstein polaron interacting with an arbitrary phonon subsystem. We assume only one electron Wannier state |r) per 
unit cell, leaving the multi-orbital case for future work. In the spirit of the Holstein model, we further assume that 
the electron interacts only with the ion degrees of freedom of the unit cell it currently occupies, i.e. 



(44) 



where r(r) is the electron coordinate at imaginary time r. Note that the local force constants Kg are different for 
different displacements. The locality of the interaction eliminates the summations over n and n' in equations (|34() - (|35|l . 
As a result, the Holstein polaron action becomes 



^o[r(r)] = ^E 



N 



kA 



kA 



Iq Jq sinh?iti>kA 



r'\) 



E 



-M 



"sA^s'A^ 



-r(r')] 



(45) 



AAAr[r(r)] 



kA 



E 



drdr 2 ^ 

sinh ?iWkA/3 



'■sA^s 



ik[r(r)- r(r')+Ar] _ ik[r(T)-r( 



(46) 



The three cases considered in this section are important on their own right, but they also serve as examples of how 
particular models are derived from the general polaron action. Many more models can be developed. In particular, 
one can envision combinations of the above cases: Holstein interaction with one dispersive phonon branch, Holstein 
interaction with uncoupled three-dimensional oscillators, and so on. 



VI. DISCUSSION 



In this paper, the polaron action for the general dispersive phonon system has been derived. It is important to list 
the conditions under which these results have been obtained. First of all, the crystal lattice is treated in the quadratic 
approximation, i.e. the phonons do not interact with each other. Secondly, the electron-phonon interaction is limited 
to the electron density - ion displacement type. There are other important classes of electron-phonon interactions - 
most notably the Jahn- Teller and Su-Schrieffer-Heeger ones - which require separate analysis. Finally, the complexity 
of the electron kinetic energy has not been fully addressed. This generalization is yet to be addressed in the polaron 
literature. 

The motivation for the present work was the extension of the capabilities of the path-integral polaron Monte Carlo^. 
However, the practical usefulness of the developed formalism entirely depends on the efficient calculation of the action. 
Consider the general formulas 1)34(1 and H35|) . Not only they involve double integration over the imaginary time, but 
also summation over all phonon momenta and branches, and a double sum over the unit cells and ion degrees of 
freedom. Thus direct calculation of the ful action at every Monte Carlo update may be prohibitively expensive. 
Development of the practical ways to apply the present formalism to specific polaron models will be the subject of 
future investigations. 
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